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1.1 Introduction 



Information processing by a network of dynamical elements is a delicate matter: 
Avalanches of activity can die out if the network is not connected enough or if the 
elements are not sensitive enough; on the other hand, activity avalanches can grow 
and spread over the entire network and override information processing as observed 
in epilepsy. 

Therefore, it has long been argued that neural networks have to establish and 
maintain a certain interme diate level of activity in order to keep away from the 
regin i es of chaos and silence ( Langton , 1990 ; Herz and Hopfield , 19951 : Bak and Chialvol . 



200 ll : iBornholdt and Rohll . l2003l ) . Similar ideas were also formulated in the context 



of genetic networks where Kauffman postulated that information processing in these 
evolved biochemical networks would be optimal near the "edge of chaos" , or cri ti- 
cality, of the dynamical percolation transition of such networks ( Kauffman . 19931 ). 
In the wake of self-organized criticality (SOC), i t was asked if al so neural sys- 



tems were self-organized to some form of criticality ( Bak et al. . 19881 ). In addition 



actual observations of neural oscillations within the hu r nan b rain were related to 
a possible SOC phenomenon (jLinkenkaer-Hansen et al.l . l200ll ). An early example 



of a SOC model that had be en adapted to be applicable to neural networks is the 
model bv lEurich etlH (|2002l l. They show that their model of the random neighbor 
Olami-Feder-Christensen universality class exhibits (subject to one critical coupling 
parameter) distributions of avalanche sizes and durations which they postulate could 
also occur in neural systems. 

Another early exa i nple of a model for self-organized critical neural networks 
(jBornholdt and Rohll . l200ll . l2003l ) drew on an alternative app r oach to self-organized 
criticality based in dynamical networks ( Bornholdt and Rohll . 200d ). Here networks 
are able to self-regulate towards and maintain a critical system state, via simple local 
rewiring rules which are plausible in the biological context. 



In "Criticality in Neural Systems", Niebur E, Plenz D, Schuster HG (eds.) 2013 (in press). 
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1 Self-organized criticality in neural network models 



After these first model approaches, indeed strong evidence for criticality in neural 
systems has bee n found in terms o f spat io-temporal activity avalanches, first in the 
seminal work of Beggs and Plenz ( 20031 ). Much further experimental evidence has 
been found since, which we will briefly review below. These experimental findings 
sparked intense research on dynamical models for criticality and avalanche dynamics 
in neural networks, which we also give a brief overview below. While most models 
emphasized biological and neurophysiological detail, our path here is different: 

The purpose of this review is to pick up the thread o f the early self-organized 
critical neural network model dBornholdt and Rohl booi ) and test its applicability 
in the light of experimental data. We would like to keep the simplicity of the first 
spin model in the light of statistical ph ysics, while lifting the drawback of a spin 
formulation w.r.t. the biological system ( Rvbarsch and Bornholdt . 2012al ). We will 
study an improved model and show that it adapts to criticality exhibiting avalanche 
statistic s that compare well with experim ental data without the need for parameter 
tuning ( Rvbarsch and Bornholdt , 2012bl ). 



1.2 Avalanche dynamics in neuronal systems 



1.2.1 Experimental results 

Let us first briefly review the experimental studies on neuronal avalanche dynamics. 
In 2003, Beggs and Plenz publis hed their findings abou t a novel mode of activ- 
ity in neocortical neuron circuits ( Beggs and Plenz . 20031 ). During in-vitro experi- 
ments with cortex slice cultures of the rat where neuronal activity in terms of local 
field potentials was analyzed via a 8x8 multi-electrode array, they found evidence 
of spontaneous bursts and avalanche-like propagation of activity followed by silent 
periods of various lengths. The observed power-law distribution of event sizes in- 
dicates that the neuronal network is maintained in a critical state. In addition, 
they found that the precise spatio-temporal patterns of the aval anches are stable 
over many hours and also robust against external perturbations (jBeggs and Plenj . 
2OO4I I. They conclude that these neuronal avalanches might play a central role for 
brain functions like information storage and processing. Also during developmental 
stages of in-vitro cortex slice cultures from newborn rats, neuronal avalanches were 
found, indic ating a homeostatic evolu tion of the cortical network during postnatal 



maturation ( Stewart and Plenz . 200?! ). Moreover, also cultures of dissociated neu- 



rons were found to exhibit this type of spontaneuous activity bu rsts in different 



kinds of networks, like rat hippocampal neurons and le ech ganglia (IMazzoni et al. 



20071 '). as well as dissociated neurons from rat embryos (jPasquale et al.l . |2008| . 



Aside from these in-vitro experiments, extensive studies in-vivo have since been 
conducted. The emergence of spontaneous neu ronal avalanches has been shown in 
anaesthesized rats during cortical development (|Gireesh and Plena. |2008 | ) as well as 



in aw ake rhesus monkeys during ongoing cortical synchronization (jPetermann et al. 
2OO9I V 
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1.2 Avalanche dynamics in neuronal systems 



The biological relevance of the avalanche-like propagation of activity in conjunc- 
tion with a critical state of the neuronal network has been emphasized in several 
works rece ntly. Such network activity has proven to be optinial for maximum dynarn - 
ical range ( Kinouchi and Copelh . 2006 : Shew et al. . 20091 : Larremore et al. . 201 ll ). 



maximal information capacity and transmission capabilit y (IShew et a,l 



well as for a maximal variability of phase synchronization (|Yang et al. 



201lh as 
2012l l. 



1.2.2 Existing models 



The experimental studies with their rich phenomenology of spatio-temporal patterns 
sparked a large number of theoretical studies and models for criticality and self- 
organization in neural networks, ranging from very simple toy models to detailed 
representations of biological functions. Most of them try to capture self-organized 
behavior with emerging avalanche activity patterns, with scaling properties similar 

to the experimental power-law event size or duration distribution s. 

As briefly mentioned above, early works as Bornholdt and Rohll ( 2000l ) and Bornholdt and Rohl 
(|200ll . I2OO3I ) focus on simple mechanisms for self-organized critic al dynamics in spin 
netw orks, which also have been discussed in a wider context (^Gross and Blasiua . 

These models represent an approach aiming at utmost simplicity of the 
model, quite similar to the universality viewpoint of statistical mechanics, rather 
than faithful representations of neurobiological and biochemical detail. Nevertheless 
they are able to self-regulate towards and maintain a critical system state, manifested 
in features as a certain limit cycle scaling behavior, via simple local rewiring rules 
which are still plausible in the biological context. We will have a closer look at these 
models in the following section \1.'6\ because they provide some of the groundwork 
for current models. 

Regarding neuronal avalanches, a 2002 work by Eurich et al. investigates net- 
works of globally coupled threshold elements which are related to integrate-and-fire 
neurons. They present a model which, after proper parameter tuning, exhibits 
avalanche-like dynamics with distinctive di stributions of avalanche sizes and dura- 
tions as expected at a critical system state ( Eurich et al. . 20021 ) . 



It is notable that these models came up even be fore experimental eviden ce was 
found for the existence of neuronal avalanches by Beggs and Plenz ( 20031 ) . Un- 
derstandably, extensive studies have been done on avalanche models following their 
discovery. Again, most models have their mechanisms of self-organization motivated 
by brain plasticity. 

A 2006 model proposed by de Arcangelis et al. consists of a model electrical 
network on a square lattice, where threshold firing dynamics, neuron refractory 
inhibition and activity-dependent plasticity of the synaptic couplings, represented 
by the conductance of the electrical links, serve as the basis for self-organization. 
Neuron states are characterized by electrical potentials, which may be emitted as 
action potentials to neighboring nodes once a certain threshold has been reached. 
With these incoming currents, the neighbor sites can eventually also reach their ac- 
tivation threshold and thus activity avalanches can propagate through the network. 
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1 Self-organized criticality in neural network models 



Avalanches are triggered by external currents to specific input sites. Following the 
activation of a node, the coupling conductances are increased by a small value for 
each link which has carried a firing current. On the other hand, after complet- 
ing a whole avalanche, all couplings in the network have their conductance reduced 
by the average of the increase which has taken place before during the avalanche 
propagation. This way, those couplings which carry many signals, will effectively 
get stronger connections, while the rather inactive connections will be reduced and 
subsequently pruned from the network. I ndeed, the model evolves to a critical state 
with power-law scaling of avalanche sizes ( de Arcangelis et al. . 20061 ) . In a following 



work, the same behavior of such an adaptive model could also be observed on a 
class of scale- free networks (namely Apollonian networks), which is argued to be 
more appropria t e as a n approach to real neuronal networks than a lattice would be 
Pellegrini et"ID (|2007l V 



A related model (in terms of insertion of li nks or facilitation o f weig hts where sig- 
nals have been passed) has been proposed bv lMeisel and GrossI ((20091). The authors 
focus on the interplay between activity-driven vs. spike-time-dependent plasticity 
and link their model to a phase transition in synchronization of the network dynam- 
ics. The temporal sequence of node activations serves as the key criterion for the 
topological updates. While they do not specifically discuss avalanche-like activity 
patterns, one observes power-law distributed quantities like correlation measures or 
synaptic weights in the self-organized states which point to dynamical criticality. 

While the last three models mentioned above are set up to strengthen those cou- 
plings ov er which signa l s have been transmitted, a kind of opposite approach was pro- 
posed by Levina et al. ( 20071 ). In their model, synaptic depression after propagation 
of activity over a link - biologically motivated by the depletion of neuro-transmitter 
ressources in the synapse - is the key mechanism which drives their model to a 
critical behavior. The network consists of fully connected integrate-and-fire neurons 
whose states are described by a membrane potential. This potential is increased 
by incoming signals from other sites or by random external input, and may cause 
the site to emit signals when the activation threshold is exceeded. Following such 
a firing event, the membrane potential is reduced by the threshold value. Again, a 
single neuron starting to fire after external input, may set off an avalanche by pass- 
ing its potential to other sites, which in turn can exceed their activation threshold, 
and so on. The couplings in this model take non-discrete values, directly related to 
the biologically relevant amount of neuro-transmitter available in each synapse. In 
short, whenever a signal is passed by a synapse, its value will be decreased. The 
coupling strength is on the other hand recovering slowly towards the maximum value 
in periods when no firing events occur. The authors also extend the model to con- 
sider leaky integrate-and-fire neurons, and find a robust self-organization towards a 
critical state, ag a,in with the charact eristic power-law distribution of avalanche sizes. 
In a later work ( Levina et al. . 20091 ). the authors further investigate the nature of 
the self-organization process in their model and discuss the combination of first- and 
second-order phase transitions with a self-organized critical phase. 

Meanwhile, field-theoretic approaches to fluctuation effects helped to shed light 
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1.2 Avalanche dynamics in neuronal systems 



on th e universality classes to expect in critical neural networks (jBuice and Cowanl . 
2007! ) and the presence of SOC in non-conservative network models of leaky neurons 
were linked to the exi stence of alternating states of high vs. low activity, so-called 
up- and down-states ( Millman et al. . 2O10l ). It has been shown that anti-Hebbian 
evolution is generally capable of creating a dynamically critical network when the 
anti-Hebbian rule affects only symmetric components of the connectivity matrix. 
The anti-symmetric component remains as an indepen dent degree of freedom and 
could be useful in learning tasks (jMagnasco et al.l . l2009l ^. Another model highlights 
the importance of synaptic plasticity for a phase transition in gen eral and relates 



the e xistence of a broad critical regime to a hierarchical modularity ([Rubinov et al. 



201 ih . The biological plausibility of activity-dependent synaptic plasti city for adap- 



tive s elf-organized critical networks has further been stressed recently (jProste et al. 



20l3). Also, robustn ess of critical bra i n dyn amics to complex network topologies 



has been emphasized ( Larremore et al. . 20121 ) 



The relevance of the critic a,l state in neuronal networks for a brain function 
as learning was underlined by de Arcangelis and Herrmann ( 2010l ). where the au- 
thors find that the perfomance in learning logical rules as well as time to learn are 
strongly connected to the strength of plastic adaptation in their model, which at 
the s ame time is able to reproduc e criti cal avalanche activity. In a most recent 
work (jde Arcangelis and Herrmannl . |2012| ) , the same authors present a new variant 
of their electrical network model. They again use facilitation of active synapses as 
their primary driving force in the self-organization, but now focus more on activity 
correlation between those nodes which are actually active in consecutive time steps. 
Here, they investigate the emergence of critical avalanche events on a variety of dif- 
ferent network types, as for example regular lattices, scale-free networks, small-world 
networks or fully connected networks. 



Also most recently, iLombardi et al.l (|2012l ) investigate the temporal organization 
of neuronal avalanches in real cortical networks. The authors find evidence that 
the specific waiting time distribution between avalanches is a consequence of the 
above-mentioned up- and down-states, which in turn is closely linked to a balance 
of excitation and inhibition in a critical network. 



While the proposed organization mechanisms strongly differ between the indi- 
vidual models (see Figure 11.11 for a cartoon representation of the various mecha- 
nisms), the resulting evolved networks tend to be part of only a few fundamental 
universality classes, exhibiting e.g. avalanche statistics in a similar way as in the 
experimental data, as power-law distributions at an exponent of —3/2. With the 
recent, more detailed models in mind, we are especially interested in the underlying 
universality of self-organization, also across other fields. Considering the enormous 
i nterest in neuronal self-organization, we here come bac k to our older spin models 
(iBornholdt and Rohlj I2000I : iBornholdt and Rohl booi ) and develop a new basic 
mechanism in the light of better biological plausibility of these models. 
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1 Self-organized criticality in neural network models 
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Figure 1.1: Schematic illustration of some of the different approaches to self-organization 
in neural network models. Rows 1-2: Links are cither added (denoted by 
-|-1; green link) or removed (denoted by -1; red link) as a function of node 
activity or correlation between nodes. Rows 3-4: Here, activity or inactivity 
of a node affects all outgoings links (thin lines). All weights of the outgoing 
links from a node are decreased (red) or increased (green) as a function of 
node activity. Row 5: Links are created and facilitated when nodes become 
active in the correct temporal sequence. Links directed against the sequence of 
activation are deleted. Row 6: Positive correlation in the activity between two 
nodes selectively increases the corresponding link, whereas there is non-selective 
weight decrease for links between uncorrelated or inactive nodes. 
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1.3 Simple models for self-organized critical adaptive neural networks 



1.3 Simple models for self-organized critical adaptive neural 
networks 

1.3.1 A first approach: node activity locally regulates connectivity 



In th e very minimal model of a random threshold network (jBornholdt and Rohlf 



a simple local mechanism for topology evolution based on node activity has 



been defined, which is capable of driving the network towards a critical connectivity 
of Kc = 2. Consider a network composed of randomly and asymmetrically con- 
nected spins (fjj = ±1), which are updated synchronously in discrete time steps via 
a threshold function of the inputs they receive: 

+ = sgn(/i(t)) (1.1) 

using 

sgn(x) = I ^ ^ Q (1.2) 

and 

N 

Mt) = Y,Wit) + h- (1.3) 

i=i 

where the link weights have discrete values Cij = ±1 (or Cij = if node i does not 
receive input from node j). In the minimal model, activation thresholds are set to 
h = for all nodes. A network run is started with random initial configuration and 
iterated until either a fixed point attractor or a limit cycle attractor is reached. The 
attractor of a network is where its dynamics ends up after a while, which is either 
a fixed point of the dynamics (all nodes reach a constant state) or a limit cycle of 
the whole network dynamics. A limit cycle in these discrete dynamical models is a 
cyclic sequence of a finite number of activation patterns. 

For the topological evolution, a random node i is selected and its activity during 
the attractor period of the network is analyzed. The network is observed until such 
an attractor is reached; and afterwards, activity of the single node during that period 
is measured. In short, if node i changes its state (jj at least once during the attractor, 
a random one of the existing non-zero in-links Cij to that node is removed. If, vice 
versa, C7j remains constant throughout the attractor period, a new non-zero in-link 
Cij from a random node j is added to the network. 

In one specific among several possible realizations of an adaptation algorithm, the 
average activity A{i) of node i over an attractor period from Ti to T2 is defined as 

1 

A{i) = —— <y^{t). (1.4) 

Topological evolution is now imposed in the following way: 

1. A random network with average connectivity Kini is created and node states 
are set to random values. 
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1 Self-organized criticality in neural network models 



2. Parallel updates according to (jl.ip are executed until a previous state reap- 
pears (i.e. until a dynamical attractor is reached). 

3. Calculate A{i) for a randomly selected node i. If \A{i)\ = 1, node i receives a 
new in-link of random weight Cij = ±1 from a random other node j. Otherwise 
(i.e. if the state of node i changes during the attractor period), one of the 
existing in-links is set to zero. 

4. Optional: a random non-zero link in the network has its sign reversed. 




Figure 1.2: Evolution of the average connectivity K^v with an activity-driven rewiring, 
shown for two different initial connectivities. Independent of the initial con- 
ditions chosen at random, the networks evolve to an average connectivity 
Kev ~ 2.55 ± 0.04. Time t is in simulation steps. 

A typical run of this algorithm will result in a connectivity evolution as shown in 
Figure [L2] for a network of = 1024 nodes. Independent from the initial connectiv- 
ity Kini, the system evolves towards a statistically stationary state with an average 
evolved connectivity of Kev{N = 1024) = 2.55 it 0.04. With increasing system size 
A, Kev converges towards Kc = 2 for the large system limit N ^ oo with a scaling 
relationship 

i^e.(A) =2 + ciV-^ (1.5) 

with c = 12.4 ± 0.5 and 6 = 0.47 ± 0.01 (compare Figure O]). 

The underlying principle which facilitates self-organization in this model is based 
on the activity A{i) of a node being closely connected to the frozen component of the 
network ~ the fraction of nodes which do not change their state along the attractor 
- which also undergoes a transition from a large to a vanishing frozen component at 
the critical connectivity. At low network connectivity, a large fraction of nodes will 
likely be frozen, and thus receive new in-links once they are selected for rewiring. 
On the other hand, at high connectivity, nodes will often change their state and in 
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1.3 Simple models for self-organized critical adaptive neural networks 
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Figure 1.3: A finite size scaling of tlie evolved average connectivity Kev as a function of 
network size N reveals that for large TV, the mean K^v converge towards the 
critical connectivity Kc = 2. For systems with N < 256, the average was taken 
over 4 x 10^ time steps, for N = 512 and N = 1024 over 5 x 10^ and 2.5 x 10^ 
time steps respectively. 



turn lose in-links in the rewiring process. Figure 11.41 shows the above mentioned 
transition as a function of connectivity for two different network sizes. With a finite 
size scaling of the transition connectivities at the respective network sizes, it can be 
shown that for large N oo, the transition occurs at the critical value of Kc = 2. 

1.3.2 Correlation as a criterion for rewiring: self-organization on a spin 
lattice neural network model 



The models described in the following section and originally proposed bv lBornholdt and Rohl 



capture self-organized critical behavior on a two-dimensional spin lattice 



via a simple correlation-based rewiring method. The m otivation behind the new ap - 
proach was to transfer the idea of self-organization from Bornholdt and Rohlf (2000) 



to neural networks, thus creating a first self-or ganized critical neura l netw ork model. 

In contrast to the activity-regulated model jBornholdt and Eohli [ioool ') discussed 
above, now: 

• the topology is constrained to a squared lattice, 

• thermal noise is added to the system, 

• link weights take non-discrete values, and 

• activation thresholds may vary from 0. 

The model consists of a randomly and asymetrically connected threshold network 
of N spins (cTj = ±1), where links can only be established among the eight local 
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1 Self-organized criticality in neural network models 
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Figure 1.4: The frozen component C(K,N) of random threshold networks, as a function of 
the networks' average connectivities K, shown for two different network sizes 
N. Sigmoid function fits can be used for a finite size scaling of the transition 
connectivity from active to frozen networks. The results indicate that this 
transi tion takes place at Kc = 2 fo r large N, details can be found in the original 
work (jBornholdt and Rohlilioooh . 



neighbors of any lattice site. The link weights Wij can be activating or inhibiting 
and are chosen randomly from a uniform distribution Wij € {— 1,+1}. The average 
connectivity K denotes the average number of non-zero incoming weights. All nodes 
are updated synchronously via a simple threshold function of their input signals from 
the previous time step: 

Prob[ai(t + l) = +1] = g^im) 

Prob[cTi(t + l) = -1] = l-g^{fi{t)) (1.6) 

with 

AT 

fi{t) = ^w,jaj{t) + Q, (1.7) 

and 

= l+exp('2^/.W) (l'^' 

where /? denotes the inverse temperature and 0j is the activation threshold of node 
i. Thresholds are chosen as @i = —0.1 + 7 where 7 is a small random Gaussian noise 
of width e. The model per se exhibits a percolation transition under variation of 
if or 0, changing between a phase of ordered dynamics, with short transients and 
limit cycle attractors, and a chaotic phase with cycle lengths scaling exponentially 
with system size. 

On a larger time scale, the network topology is now changed by a slow local 
rewiring mechanism according to the following principle: if the dynamics of two 
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1.3 Simple models for self-organized critical adaptive neural networks 

neighboring nodes is highly correlated or anti-correlated, they get a new link between 
them. Otherwise, if their activity shows low correlation, any present link between 
them is removed, which is reminiscient of the Hebbian learning rule. In this model, 
the correlation Cjj(T) of nodes i,j over a time interval r is defined as 

C^J{r) = -—Y,^^{t)<y,{t). (1.9) 

t=to 

The full model is constructed as follows: 

1. Start with a randomly connected lattice of average connectivity Kini, random 
initial node configuration, and random individual activation thresholds 0j. 

2. Perform synchronous updates of all nodes for r time steps (Here, the choice 
of T is not linked to any attractor period measurement, but should be chosen 
large enough to ensure a separation of time scales between network dynamics 
and topology changes.). 

3. Choose a random node i and random neighbor j and calculate Cij{T/2) over 
the last r/2 time steps (a first r/2 time steps are disregarded to exclude a 
possible transient period following the previous topology change.). 

4. If |Cij(r/2)| is larger than a given threshold a, a new link from node j to i is 
inserted with random weight Wij S {—1, +1}- If else |Cij(r/2)| < a the weight 
Wij is set to 0. 

5. Use the current network configuration as new initial state and continue with 
step 2. 




0.01 



0.001 




Figure 1.5: Left: Evolution of the average connectivity with a correlation-based rewiring, 
shown for two different initial connectivities. Again, connectivity evolves to a 
specific average depending on network size, but independent from the initial 
configuration. Right: Finite size scaling of the evolved average connectivity. 
The best fit is obtained for a large system limit of = 2.24 ±0.03. Averages 
are taken over 4 x 10^ time steps. 
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1 Self-organized criticality in neural network models 

Independent from the initial average connectivity Kmi, one finds a slow conver- 
gence of K towards a specific mean evolved connectivity K^y , which is characteristic 
for the respective network size N (Figure I1.5P and shows a finite size scaling accord- 
ing to 

Kev{N) = aN~^ + b (1.10) 

with a = 1.2 ± 0.4, 5 = 0.86 ± 0.07 and b = 2.24 ± 0.03, where b can be interpreted 
as the evolved average connectivity for the large system limit N ^ oo: 

iCg"^ = 2.24 ± 0.03 (1.11) 




5 10 15 20 25 30 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



/9 mr]\ 

Figure 1.6: Left: Evolved average connectivity K^v as a function of the inverse temperature 
(3. The behavior is robust over a wide range of /?. Each point is averaged over 
10^ time steps in a network of size N = 64. Right: Histogram of the average 
correlation | (r) | for a network evolving in time with = 64 and f3 =- 10. As 
very low and very high correlations dominate, the exact choice of the correlation 
threshold during the rewiring process is of minor importance. 

In addition, it is shown that the proposed adaptation mechanism works robustly 
towards a wide range of thermal noise /?, and also the specific choice of the corre- 
lation threshold a for rewiring only plays a minor role regarding the evolved Kf.^ 
(Figure [L6]). 

Having a closer look at a finite size scaling of the evolved average attractor length, 
one finds a scaling behavior close to criticality. While attractor lengths normally 
scale exponentially with system size in the supercritical, chaotic regime and sublin- 
early in the subcritical, ordered phase, this model exhibits relatively short attractor 
cycles also for large evolved networks in the critical regime (Figure [1.7|) . 

1.3.3 Simplicity vs. biological plausibility - and possible improvements 
Transition from spins to Boolean node states 

In the above sections, we have seen that already basic toy models, neglecting a lot 
of details, can reproduce some of the observations made in real neuronal systems. 
We now want to move these models a step closer towards biologigal plausibility and 
at the same time construct an even simpler model. 
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1.3 Simple models for self-organized critical adaptive neural networks 




Figure 1.7: Finite size scaling of the evolved average attractor period A{N) for networks of 
different sizes N. (b). Also shown for comparison is the corresponding scaling 
of the attractor lengths of a supercritical random network (a) with K ~ 3.8 
and a subcritical one (c) with K = 1.5. 

One major shortcoming of both models discussed above is the fact that they are 
constructed as spin models. In some circumstances, however, when faithful repre- 
sentation of certain biological details is important, the exact definition matters. In 
the spin version of a neural network model, for example, a node with negative spin 
state (jj = —1 will transmit non-zero signals through its outgoing weights Cij, de- 
spite representing an inactive (!) biological node. In the model, such signals arrive 
at target nodes i, e.g., as a sum of incoming signals fi = '^jLiCijaj. However, 
biological nodes, as genes or neurons, usually do not transmit signals when inac- 
tive. Also in other contexts like biochemical network models, each node represents 
whether a specific chemical component is present {a = 1) or absent (a = 0). Thus 
the network itself is mostly in a state of being partially absent as, e.g., in a protein 
network where for every absent protein all of its outgoing links are absent as well. In 
the spin state convention, this fact is not faithfully represented. A far more natural 
choice would be to construct a model based on Boolean state nodes, where nodes 
can be truly "off' (ai = 0) or "on" (ui = 1) - which is precisely what we are going 
to do in the following sections. 

Another example for an inaccurate detail is the common practice to use the stan- 
dard convention of the Heaviside step function as an activation function in discrete 
dynamical networks (or the sign function in the spin model context). The conven- 
tion 0(0) = 1 is not a careful representation of biological circumstances. Both, for 
genes and neurons, a silent input frequently maps to a silent output. Therefore, we 
use a redefined threshold function defined as 

e«W = {o': x<a 

Most importantly, in our context here, the choice of Boolean node states and 
the redefined threshold function are vital when we discuss mechanisms of self- 
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1 Self-organized criticality in neural network models 



organization . For instance, the correla tion-based approach presented in the older 
spin model (jBornholdt and R,5hi booi ) in section 0:2] explicitly measures contr 



1- 

butions by pairs of nodes which could be constantly off (cjj/j = —1) and still treat 
them as highly correlated (because (—1) • (—1) = +1) even though there is no ac- 
tivity involved at all. In the later section 11.3. 4^ we will therefore present a new 
approach for a network of Boolean state nodes, which does not rely on non-activity 
correlations anymore. 

Model definitions 

Let us consider randomly wired threshold networks of N nodes G {0, 1}. At each 
discrete time step, all nodes are updated in parallel according to 

ai{t + 1) = eoifiit)) (1.13) 

using the input function 

N 

fi{t) = Y,w{t) + ei. (1.14) 

In particular we choose 0o(O) := for plausibility reasons (zero input signal will 
produce zero output). While the weights take discrete values Cij = ±1 with equal 
probability for connected nodes, we select the thresholds 6i = for the following 
discussion. For any node i, the number of incoming links Cij 7^ is called the 
in-degree ki of that specific node. K denotes the average connectivity of the whole 
network. With randomly placed links, the probability for each node to actually have 
incoming links follows a Poissonian distribution: 

p{ki = k) = ^-exp{-K). (1.15) 

Exploration of critical properties - Activity-dependent criticality 

To analytically derive the critical connectivity of this type of network model, we first 
study damage spreading on a local basis and calculate the probability Ps{k) for a 
single node to propagate a small perturbation, i.e. to change its output from to 1 
or vice versa after changing a single input state. The calculati on can be done closely 



follow ing the derivation for spin-type threshold networks by iRohlf and Bornholdt 



(j2002l ). but one has to account for the possible occurrence of '0' input signals also via 
non-zero links. The combinatorial approach yields a result that directly corresponds 
to the spin-type network calculation viap^°°'(fc) = p|^™(2/c). However, this approach 
does not hold true for our Boolean model in combination with the defined Theta 
function 0o(O) := as it assumes a statistically equal distribution of all possible 
input configurations for a single node. In the Boolean model, this would involve 
an average node activity of 6 = 0.5 over the whole network (where b denotes the 
average fraction of nodes which are active, i.e. = 1). Instead we find (Fig. II. 8p 
that the average activity on the network is significantly below 0.5. At X = 4 (which 
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Figure 1.8: Average node activity b as function of connectivity K measured on attractors 
of 10000 sample networks each, 200 nodes. 



will turn out to be already far in the supercritical regime), less than 30 percent of 
all nodes are active on average. Around K ^ 2 (where we usually expect the critical 
connectivity for such networks), the average activity is in fact below 10 percent. 
Thus, random input configurations of nodes in this network will more likely consist 
of a higher number of '0' signal contributions than of ±1 inputs. 

Therefore, when counting input configurations for the combinatorial derivation 
of Ps{k), we need to weight all relevant configurations according to their realiza- 
tion probabili ty as given by the average activ ity b - the detailed derivation can 
be found in ( Rvbarsch and Bornholdt . 2012al ). With the average probability of 



damage spreading, we can further compute the branching parameter or sensitiv- 



ity A = K ■ (p,c,)(K) and apply th e annealed approximation (jPerrida and Pomeau 



198fil : iBornholdt and Rohli I2OO0I ) to obtain the critical connectivity Kc by solving 

Xc = K,-{ps){K,) = l. (1.16) 

However, Kc now depends on the average network activity, which in turn is a 
function of the average connectivity K itself as shown in Fig. 11.81 Prom the combined 
plot in Pig. 11.91 we find that both curves intersect at a point where the network 
dynamics - due to the current connectivity K - exhibit an average activity which 
in turn yields a critical connectivity Kc that exactly matches the given connectivity. 
This intersection thus corresponds to the critical connectivity of the present network 
model. 

However, the average activity still varies with different network sizes, which is ob- 
vious from Figure [L9l Therefore, also the critical connectivity is a function of A^. For 
an analytic approach to the infinite size limit, it is possible to calculate the average 
network activity at the next time step bt+i in an iterative way from the momentary 



averag e activity 6f Again, the details can be found in (jRvbarsch and Bornholdt 
2012a|). We can afterwards distinguish between the different dynamical regimes by 
solving (fet+i) = bt{K) for the critical line. The solid line in Figure fTTU] depicts the 
evolved activity in the long time limit. We find that for infinite system size, the 
critical connectivity is at 

Kc{N ^ 00) = 2.000 ± 0.001 



15 



1 Self-organized criticality in neural network models 




Figure 1.9: Average activity b on attractors of different network sizes (right to left: N = 
50, 200, 800, ensemble averages were taken over 10000 networks each). Squares 
indicate activity on an infinite system determined by finite size scaling, which is 
in good agreement with the analytic result (solid line). The dashed line shows 
the analytic result for Kc{b) from eq. (|1.16|) . The intersections represent the 
value of Kc for the given network size. 



while up to this value all network activity vanishes in the long time limit {b^o = 0). 
For any average connectivity K > 2, a, certain fraction of nodes remains active. 
In finite size systems, both network activity evolution and damage propagation 
probabilities are subject to finite size effects, thus increasing Kc to a higher value. 
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Figure 1.10: Average attractor length at different network sizes. Ensemble averages were 
taken over 10000 networks each at (A) K = 2.4, (B) K = 2.0, (C) K = 1.6. 
Inset figure shows the scaling behavior of the corresponding transient lengths. 

Finally, let us have a closer look on the average length of attractor cycles and 
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transients. As shown in Fig. I1.10| the behavior is strongly dependent of the dy- 
namical regime of the network. As expected and in accordance with early works 
on random threshold netw orks ( Kiirten . 19881 ) as well as random Boolean networks 



( Bastolla and Parisi . 19961 ) . we find an exponential increase of the average attractor 



lengths with network size in the chaotic regime {K > Kc), whereas we can observe 
a power-law increase in the frozen phase {K < Kc)- We find similar behavior for 
the scaling of transient lengths (inset of Figure ll.lOp . 

Extension of the model: thermal noise 

As it is clear from eq. ()1.13p . nodes in our model will only switch to active state if 
they get a positive input from somewhere. Thus, to get activity into the system, 
we could either define certain nodes to get an additional external input, but this 
would at the same time create two different kinds of nodes, those with and those 
without external activity input, which would in turn diminish the simplicity of our 
model. That is why we will alternatively use thermal noise to create activity, using 
a Glauber update of the no des in the same way as it was discussed in the spin model 
()Bornholdt and Rohilionii ) in section [TMl with one slight modification. We define 



a probability for the activation of a node, which is a sigmoid function of the actual 
input sum, but also leaving room for a spontaneous activation related to the inverse 
temperature (3 of the system. 

Prob[a,(t + l) = l] = g^im) 

Prob[a,(t + l) = 0] = 1-gpim) (1.17) 

with 

N 

i=i 

and ^ 

1 -hexp(-2/^(/j(t) -0.5)) 



You will note the similarity to the older spin model (|Bornholdt and Rohj \200i ) 



but be aware that in eq. ()1.19p we shift the input sum fi by —0.5, and we will explain 
now why this is necessary. Remember that we use binary coupling weights Cij = ±1 
for existing links in our model. The input sum fi to any node will therefore be an 
integer value. If we would not shift the input sum, a node with an input fi = 
(which should always be inactive in the deterministic model without thermal noise), 
would, after the introduction of the Glauber update rule with non-zero temperature, 
always have a probability of Prob[(Tj(t -|- 1) = -|-1] = 0.5 to be activated, regardless 
of the actual inverse temperature /3. Figure [1.111 illustrates this problem. A simple 
shift of —0.5 will now give us the desired behavior: Nodes with input fi = will stay 
off in most cases, with a slight probability to turn active depending on temperature, 
and on the other hand nodes with activating input fi = +l will be on in most cases, 
with slight probability to remain inactive. With this modification of the original 
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Figure 1.11: With integer input sum values, nodes without input would turn active with 
probability 0.5 regardless of temperature. To prevent this, we shift the input 
sum by —0.5 such that the transition between off and on state happens exactly 
in the middle between input sums resp. 1. 



basic model ( Rvbarsch and Bornholdt . 2012al ). we can now continue to make our 



network adaptive and critical. 

1.3.4 Self-organization on the Boolean state model 

We now want to set up an adaptive network based on the model discussed above, 
which is still capable of self-regulation towards a critical state despite being simplified 
to the most minimal model possible. Particularly, we want it to 

• have a simple, yet biologically plausible rewiring rule, which only uses local 
information accessible to individual nodes 



work independently from a special topology as a lattice. 



We have a,lready pointed out the prob lems of a spin formulation of neural network 
models like Bornholdt and Rphil hmis^ . and possible ways out with the Boolean 



state model. As a major advantage of the latter, activity avalanches intrinsically 
occur in this type of network, whereas spin networks typically exhibit continuous 
fluctuations with no avalanches directly visible. However, the old correlation-based 
rewiring mechanism will no longer work when inactive nodes are now represented 
by '0' instead of '—1'. A solution will be presented below. 

A second aspect which needs to be addressed concerning the self-organization 
mechanism is topological invariance of the algorithm. The older, correlation-based 
mechanism from the spin model relies on randomly selecting neighboring sites on a 
lattice for the rewiring processes. On a lattice, the number of possible neighbors is 
strictly limited, but on a large random network near critical connectivity, there are 
far more unconnected pairs of nodes than there are connected pairs. Thus, randomly 
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selecting pairs of nodes for rewiring would introduce a strong bias towards connect- 
ing nodes which were previously unconnected. This results in a strong increase of 
connectivity and makes a self-organized adaptation towards a critical state practi- 
cally impossible. If we want to overcome the restriction of e.g. a confined lattice 
topology in order to improve biological applicability of the model, we have to adapt 
the rewiring mechanism such that this bias no longer exists. 

Of course, the new model shall inherit several important features from the older 
spin models, which already underline the applicability to biological networks: in 
particular, it must be capable of self-regulation towards a critical state despite being 
simplified to the most minimal model possible. The organization process however 
should be based on a simple, yet biologically plausible rewiring rule, which only uses 
local information accessible to individual nodes like pre- and post-synaptic activity 
and correlation of such activity. 

This section is directly based on our recent work ( Rvbarsch and Bornholdt . 2012bl ). 



Model definitions 

Consider a randomly connected threshold network of the type discussed above in 
section ri.3.31 The network consists of N nodes of Boolean states € {0, 1} which 
can be linked by asymmetric directed couplings Cjj = ±1. Node pairs which are not 
linked have their coupling set to Cij = 0; and links may exist between any two nodes, 
so there is no underlying spatial topology in this model. 

All nodes are updated synchronously in discrete time steps via a simple threshold 
function of their input signals with a little thermal noise introduced by the inverse 
temperature /?, in the same way as in the model of Bornholdt and Roehl (section 
ll.3.2p . but now with an input shift of —0.5 in the Glauber update, representing the 
new 00 function as discussed in section ll.3.3t 

Prob[a,(t + l) = l] = gfiim)) 

Prob[ai(t + 1) = 0] = l-g^{fi{t)) (1.20) 

with 

N 

m = J2ci,aj{t) + ei (1.21) 

and 

''^^'^'^^ = l + exp(-2/3(/.(t)-0.5))- ^'-''^ 

For the simplicity of our model, we first assume that all nodes have an identical 
activation threshold of 0j = 0. 

Rewiring algorithm 

The adaptation algorithm is now constructed in the following way: We start the 
network at an arbitrary initial connectivity Kini and do parallel synchronous updates 
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on all nodes according to eq. (jl.l7p . All activity in this model originates from small 
perturbations by thermal noise, leading to activity avalanches of various sizes. In 
this case we set the inverse temperature to /3 = 5. On a larger time scale, i.e. 
after r = 200 updates, a topology rewiring is introduced at a randomly chosen, 
single node. The new element in our approach is to test whether the addition or 
the removal of one random in-link at the selected node will increase the average 
dynamical correlation to all existing inputs of that node. By selecting only one 
single node for this procedure, we effectively diminish the bias of selecting mostly 
unconnected node pairs - but retain the biologically inspired idea for a Hebbian, 
correlation-based rewiring mechanism on a local level. 

Now, we have to define what is meant by dynamical correlation in this case. 
We here use the Pearson correlation coefficient to first determine the dynamical 
correlation between a node i and one of its inputs j: 

^ {aijt + l)ajit)) - {aijt + l))(cT,-(t)) ^3) 
Si ■ Sj 

where Si and Sj in the denominator depict the standard deviation of states of the 
nodes i and j, respectively. In case one or both of the nodes are frozen in their state 
(i.e. yield a standard deviation of 0), the Pearson correlation coefficient would not 
be defined, we will assume a correlation of Cij = 0. Note that we always use the 
state of node i at one time step later than node j, thereby taking into account the 
signal transmission time of one time step from one node to the next one. Again, 
as in the model of Bornholdt and E,5hll tooA ). the correlation coefficient is only 



taken over the second half of the preceding r time steps in order to avoid transient 
dynamics. Finally, we define the average input correlation C^^^ of node i as 

1 ^ 

Cr = ^T.\'^J\C^J (1.24) 

' j=o 

where ki is the in-degree of node i. The factor \cij\ ensures that correlations are 
only measured where links are present between the nodes. For nodes without any 
in-links {ki = 0) we define C'^^^ := 0. 

In detail, the adaptive rewiring is now performed in the following steps: 

1. Select a random node i and generate three clones of the current network topol- 
ogy and state: 

Network 1: This copy remains unchanged. 

Network 2: In this copy, node i will get one additional in-link from a randomly 
selected other node j which is not yet connected to i in the original copy 
(if possible, i.e. ki < N). Also the coupling weight Cij it 1 of the new link 
is chosen at random. 

Network 3: In the third copy, one of the existing in-links to node i (again if 
possible, i.e. ki > 0) will be removed. 
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2. All three copies of the network are individually run for r = 200 time steps. 

3. On all three networks, the average input correlation Cf^^ of node i to all of 
the respective input nodes in each network is determined. 

4. We accept and continue with the network which yields the highest absolute 
value of C^^^ , the other two clones are dismissed. If two or more networks 
yield the same highest average input correlation such that no explicit decision 
is possible, we simply continue with the unchanged status quo. 

5. A new random node i is selected and the algorithm starts over with step 1. 
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(2012) 





CORRELATED 



CORRELATED 



Figure 1.12: Schematic illustration of the rewiring mechanism based on average input cor- 
relation. In this example, the target node initially has three in-links. Left: If 
the addition of a fourth input increases the average input correlation C""^, a 
link will be inserted. Right: If removal of an existing in-link increases C™", 
the Unk will be deleted. 



It is worth noting that this model - in the same way as the earlier work by 
Bornholdt and R6hl (|2003l ') - uses locally available information at synapse level and 
takes into account both pre- and post -synaptic activity. This i s a fundamental 
differe nce to approache s discu ssed e.g. by de Arcangelis et al. ( 20061 ). Pellegrini et al 



(j2007l ) or iLevina et al.l (j2007l ) , where only pre-synaptic activity determines changes 
to the coupling weights. 

Note that the non-locality of running three network copies in parallel that we use 
here is not critical. A local implementation is straightforward, locally estimating 
the order parameter (here the average input correlation C^^^) as time average, with 
a sufficient time scale separation towards the adaptive ch anges in the network. A 
l ocal ve rsion of the model is studied in the latest version of lRvbarsch and Bornholdt 
(|2012bl l. 



Observations 

In the following, we will have a look at different observables during numerical simu- 
lations of network evolution in the model presented above. Key figures include the 
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average connectivity K and the branching parameter (or sensitivity) A. Both are 
closely linked to and influenced by the ratio of activating links p, which is simply 
the fraction of positive couplings Qj = +1 among all existing (non-zero) links. 
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Figure 1.13: Regardless of initial connectivity and sensitivity, the network evolves to a 
critical configuration. Left: starting with completely disconnected nodes and 
obviously subcritical "network". Right: starting with a supercritical network. 



The left part in Figure 11.131 shows an exemplary run of the topology evolution al- 
gorithm, where we start with completely isolated nodes without any links. Trivially, 
the "network" is subcritical at this stage, which can be seen from the branching 
parameter which is far below 1. As soon as rewiring takes place, the network starts 
to insert new links, obviously because these links enable the nodes to pass signals 
and subsequently act in a correlated way. With increasing connectivity, also the 
branching parameter rises, indicating that perturbations start to spread from their 
origin to other nodes. When the branching parameter approaches 1, indicating that 
the network reaches a critical state, the insertion of new links is cut back. The 
processes of insertion and depletion of links tend to be balanced against each other, 
regulating the network close to criticality. 

If we, on the other hand, start with a randomly interconnected network at a 
higher connectivity like Kmi = 4 as shown in the right-hand side of Figure II. 13^ 
we find the network in the supercritical regime (A > 1) at the beginning. When 
above the critical threshold, many nodes will show chaotic activity with very low 
average correlation to their respective inputs. The rewiring algorithm reacts in the 
appropriate way by deleting links from the network, until the branching parameter 
approaches 1. 

In both examples above, the evolution of the ratio of activating links p (which 
tends towards 1) shows, that the rewiring algorithm in general favors the insertion 
of activating links and vice versa the deletion of inhibitory couplings. This appears 
indeed very plausible if we remind ourselves that the rewiring mechanism optimizes 
the inputs of a node towards high correlation on average. Also, nodes will only switch 
to active state o"i = 1 if they get an overall positive input. As we had replaced spins 
by Boolean state nodes, this can only happen via activating links - and that is why 
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correlations mainly originate from positive couplings in our model. As a result, we 
observe the connectivity evolving towards one in-link per node, with the ratio of 
positive links also tending towards one. 
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Figure 1.14: Network evolution with activating links ratio kept at p = 0.8. 



For a richer pattern complexity, we might later want to introduce a second mech- 
anism which balances out positive via negative links, and with a first approach we 
can already test how the rewiring strategy would fit to that situation: if, after each 
rewiring step, we change the sign of single random links as necessary to obtain a 
ratio of e.g. 80% activating links (i.e. p = 0.8), keeping the large majority of present 
links unchanged, the branching parameter will again stabilize close to the critical 
transition, while the connectivity is maintained at a higher value. Figure IT . 141 shows 
that the self-organization behavior is again independent from the initial conditions. 
This result does not depend on the exact choice of the activating links ratio p; sim- 
ilar plots can easily be obtained for a large range starting ai p = 0.5, where the 
connectivity will subsequently evolve towards a value slightly below K = 2, which is 
the value we would expect as the critical connectivity for a randomly wired network 
with balanced link ra tio according to the calcu l ations made in section 11.3.31 for the 
basic network model ( Rvbarsch and Bornholdt . 2012al ). 

In addition to the branching parameter measurement, we also take a look at some 
dynamical properties of the evolved networks to further verify their criticality. As 
stated in the introduction section, we are especially interested in the resulting activ- 
ity avalanches on the networks. Figure [T. 151 (left) shows the cumulative distribution 
of avalanche sizes in an evolved sample network of = 1024 nodes. We observe a 
broad distribution of avalanche sizes and a power-law like scaling of event sizes with 
an exponent close to —1/2, corresponding to an exponent of —3/2 in the probability 
density - which is characteristic of a critical branching process. At the same time, 
this is in perfect agreement with the event size distribution observed by Beggs and 
Plenz in their in-vitro experiments. 

If we randomly activate small fractions of nodes in an otherwise silent network 
(single nodes up to 5 % of the nodes) to set off avalanches, we can also see 
(Figure fLlSl right) that the resulting transient length shows a power-law scaling with 
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Figure 1.15: Left: cumulative distribution of avalanche sizes in an evolved sample network 
of = 1024 nodes. We find a broad, near power-law distribution comparable 
to a slope of —1/2, indicative of a critical branching process and corresponding 
nicely to the experimental results of Beggs and Plenz. Right: scaling of average 
transient lengths at different network sizes, 50 evolved networks each. 



network size right up to network snapshots taken after an evolution run at the final 
average branching parameter of one. Intermediate networks taken from within an 
evolution process at a higher branching parameter instead show a super polynomial 
increase of transient lengths with system size, which is precisely what we expect. 



1.3.5 Response to external perturbations 

In recent in-vitro experiments, IPlend (|2012l l could further demonstrate that cor- 
tical networks can self-regulate in response to external influences, retaining their 
functionality while avalanche-like dynamics persist - for example after neuronal ex- 
citability has been decreased by adding an antagonist for fast glutamatergic synaptic 
transmission to the cultures. 

To reproduce such behavior in our model, we can include variations in the ac- 
tivation thresholds Qi of the individual nodes, which had been set to zero in the 
above discussions for maximum model simplicity. Assume we start our network evo- 
lution algorithm with a moderately connected, but subcritical network, where all 
nodes have a higher activation threshold of = 1. According to the update rule 
(|1.17|) . now at least two positive inputs are necessary to activate a single node. As 
the rewiring algorithm is based on propagation of thermal noise signals, the inverse 
temperature (3 needs to be selected at a lower value than before. (As a general rule, 
(3 should be selected in a range where thermal activation of nodes occur at a low 
rate, such that avalanches can be triggered, but are not dominated by noise.) 

Figure 11.161 shows that, same as above, the subcritical network starts to grow 
new links, thereby increasing the average branching parameter. Again, this process 
is stopped after the supercritical regime is reached. While the system does not 
approach to a phase transition as nicely as shown above for activation thresholds of 
zero (in fact the branching fluctuates much more around the target value of one) , the 
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Figure 1.16: Rewiring response to a sudden decrease of activation thresholds. AU 9^ were 
set from 1 to in the same time step. 

general tendency remains: the rewiring mechanism reacts properly as the network 
drifts too far off from criticality. At one time step in the center of Figure \1.16\ 
we suddenly reset all nodes to an activation threshold of @i = 0, simulating the 
addition of a stimulant. As we can expect, this immediately puts the network into a 
supercritical, chaotic state. This is reflected by the branching parameter, which now 
constantly stays above one and does not fluctuate below anymore. It is clearly visible 
that the rewiring mechanism promptly reacts and drives back the connectivity, until 
the branching parameter slowly converges towards one again. A similar behavior is 
also found if thresholds 0j are not changed all at once, but gradually during network 
evolution. 

1.4 Conclusion 

We have seen that already very minimalistic binary neural network models are ca- 
pable of self-organized critical behavior. While older models show some drawbacks 
regarding biological plausibility originating directly from their nature as spin net- 
works, we subsequently presented a possible transition to a self-organized critical, 
randomly wired network of Boolean node states with emerging dynamical patterns, 
namely activity avalanches, reminiscent of activity modes as observed in real neu- 
ronal systems. This is possible via a simple, locally realizable, rewiring mechanism 
which uses activity correlation as its regulation criterion, thus retaining the biolog- 
ically inspired rewiring basis from the spin version. While it is obvious that there 
are far more details involved in self-organization of real neuronal networks - some 
of which are reflected in other existing models - we here observe a fundamental 
organization mechanism leading to a critical system that exhibits statistical prop- 
erties pointing to a larger class of universality, regardless of the details of a specific 
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biological implementation. 
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